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We study in detail a recently proposed simple discrete model for evolution on smooth landscapes. 
An asymptotic solution of this model for long times is constructed. We find that the dynamics of the 
population are governed by correlation functions that although being formally down by powers of 
N (the population size) nonetheless control the evolution process after a very short transient. The 
long-time behavior can be found analytically since only one of these higher-order correlators (the 
, two-point function) is relevant. We compare and contrast the exact findings derived herein with a 

■ previously proposed phenomenological treatment employing mean field theory supplemented with 

a cutoff at small population density. Finally, we relate our results to the recently studied case of 
mutation on a totally flat landscape. 
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' I. INTRODUCTION 

r- 

[-»«. , In a recent paper, |lj we introduced a model to describe evolution on a smooth landscape. This model was 
motivated by the results of recent experiments |§] on the evolution of fitness in populations of RNA viruses. In these 
experiments, the population shows secular changes in birth rate, which is identified by the biologists as the logarithm 
of fitness. These experiments showed a rapid initial increase of birth rate, followed by a sharp transition to a regime 
of significantly slower, approximately linear, increase. The initial rapid rise in birth rate was interpreted by us to be a 
result of the exponentially quick dominance of the population by the most "fit" (most rapidly reproducing) members 
of the initial population. These most "fit" individuals produce many more offspring, so they quickly dominate, in 
accord with a simple picture of the working of selection in a population. The subsequent slow rise was interpreted as 
being due to the effect of mutations which, as the mutation rate is low, act on a slower time scale. 

The relative smoothness of the increase of birth rate in this second, linear, regime indicates that the "fitness" 
landscape must be fairly smooth. It is clear from many studies of evolution [^|~|| in rough, "glassy" landscapes, that 
increases in fitness in such cases are expected to be sudden, with long intervals of little or no improvement in between. 
While it is clear that the population should smoothly "climb the hill" in a smooth landscape, no attention seems 
to have been focused on the precise dynamics of such hill-climbing. In particular, what is the speed with which the 
population climbs, and on what does it depend? 

To begin to address these questions, we formulated W a very simple discrete model of evolution in a smooth 
landscape. This model exhibited both the initial rapid "collapse" of the population onto its most fit members, and 
the subsequent smooth rise in birth rate which characterized the RNA virus experiments (see Fig 1.). In addition, 
we constructed a continuum mean-field treatment which, when suitably doctored with a cutoff, appeared to give a 
qualitatively correct description of the dynamics. As we were concerned with the gross qualitative features of the 
model and its possible relevance to the experiment, we did not investigate the precise extent of the apparent agreement 
nor indicate in any detail as to where the cutoff came from. 

In this paper, we return to the analysis of our discrete model. Our analysis provides an exact description of the 
dynamics at asymptotically long times. This asymptotic dynamics indicates that the population does indeed increase 
its birth rate in a linear fashion, and provides a prediction of the rate of such increase. It also provides a description of 
the stochastic departures from this linear behavior. The most striking result of this analysis is that the rate of increase 
of birth rate (the velocity of climbing of the fitness hill) depends crucially on the population size (increasing essentially 
linearly with the population). Naive mean-field theory is however completely independent of N, the population size. 
Indeed, we show explicitly how naive mean-field theory gives rise to a finite-time divergence (or in another variant of 
the model, an exponential blowup) of the fitness and how the true dynamics cures this unphysical behavior through 
terms which are formally lower order in 1/N. These terms become relevant on relatively short time-scales and drive 
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the system to an N dependent asymptotic state. It is important to note that even the cutoff version of mean-field 
theory, which does not suffer from the blowup phenomenon, does not correctly predict this asymptotic state, due to 
the neglect of a specific two-body correlation term which dominates at long-time. This diagnosis of the failings of 
mean-field theory is the second major product of our analysis. 
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FIG. 1. Evolution of fitness of MARM clone during the transmission series of 80 passages on HeLA cells (Fig. 2b of § with 
permission from PNAS) 

II. PRELIMINARIES 

We begin our exposition by presenting the discrete model of evolution which forms the basis of our discussions. We 
also outline the naive mean- field theory and the modified cutoff version we analyzed in our previous work jjj . 

A. The Discrete Model 

The discrete model whose behavior we will primarily focus on is a slight variant of that discussed in our earlier 
work. The original model, as well as another variant, the flat landscape model, can be easily analyzed in a manner 
parallel to that with which we shall treat the version we present here. The discussion of these other variants we will 
thus postpone to a later section. 

Our model consists of a fixed population of N individuals. Although in the experiments the population varied 
between two essentially fixed values, the fixed population constraint is obviously the simpler case to study first, and 
seems to preserve the qualitative features of the more general case. The individuals are characterized by their birth 
rates, which the experimentalists associate with the logarithm of the "fitness" of the individual. Each individual gives 
birth at his characteristic rate in a standard Poisson process. The baby, due to mutation, may have a birth rate which 
is ±A that of his parent. The overall mutation probability for a baby is a constant fi. By rescaling time by A, we 
may take the change in the baby's birth rate to be ±1, so that the birth rates are restricted to be integers. We do 
not allow a baby to mutate to birth rate, so in fact the birth rates are all natural numbers. As we work at fixed 
population size N, each birth is occasioned by the unfortunate demise of some member of the group. The victim is 
chosen at random from among all members of the population (including the new-born baby) independent of their 
fitness rating. 

To get some idea of the dynamics of this model, we show in Figure 2 the results of a simulation of this model with 
a population size oi N = 50. Starting from a population of identical individuals, the average mean-birth-rate (note: 
average means over different realizations, mean denotes the mean for all individuals in a given simulation) accelerates 
quickly, followed by a gradual deacceleration to a regime of constant velocity, i.e. a linear increase of mean-birth-rate 
with time. Looking also at the variability of the population, we calculate the average sample width-squared, which 
increases rapidly from zero and then saturates at some finite value. We will see later that the transient phase is 
complicated, but the asymptotic state is always one of constant velocity. 

B. Mean-Field Theory 

It is easy to construct a mean-field theory for the simple process outlined above. Let P(n, t) be the number of 
individuals with birthrate n at time t. Then, in the absence of mutations, it is clear that 
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P = [n- n)P(n) 



(1) 



where n is the instantaneous mean value of n. The presence of mutations means that a fraction [i of the nP(n) 
newborns will leave birth-rate n, half moving up to n + 1 and half down to n — 1. Thus, 

P(n) = (n - n)P(n) - n^P(n) + i/i [(n + l)P(n + 1) - (n - l)P(n - 1)] (2) 

or, in continuum language, replacing the discrete n by x 

P(x) = (x - x)P(x) + -/z(a:P(a:))" (3) 

This equation is essentially identical to that of the mean-field theory for Diffusion-Limited Aggregation (DLA) || in 
the far-field (large- a:) regime. It is not surprising, therefore, that both the mean-field DLA |tJ and our mean-field 
Equation (Q) share the property that an initial pulse accelerates to infinite velocity in finite time. We will return to 
our mean-field Equation (|J), and in particular demonstrate explicitly the finite-time singularity, in Section VI. 

N=50 (i=0.5 
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FIG. 2. a) Mean birth-rate for population of N = 50 individuals as a function of time, with mutation rate fj, = .5, averaged 
over 200 realizations. Dotted line shows theoretical predicted rate of increase of mean-birth-rate, 
b) Width-squared from same simulation. 

For the moment, it is important to inquire as to the origin of this unphysical behavior. From a mathematical 
standpoint, the essential problem is the term xP, which is unbounded. Physically, although it is extremely unlikely 
for an individual to find itself far in front of the rest of the population, should it happen, then that individual's growth 
rate would be enormously greater than anyone else's. Thus, this very unlikely scenario makes a significant contribution 
to the ensemble average P(x) since its contribution is the product of the (tiny) probability of this configuration times 
the (huge) number of individuals at this site. In the model, however, the contribution of such an outlyer to the 
ensemble average remains negligible since the total number of individuals in any given configuration is limited to 
N. From this discussion, however, we can surmise that the velocity of propagation should be a strongly increasing 
function of N. Since mean-field theory is a theory for infinite N, mean-field theory is in some sense giving the correct 
answer, which is an infinite velocity at infinite N. While perhaps correct in this sense, such an answer is of course of 
little use. 

There are various prescriptions one might employ to remove this ugly feature of the mean-field theory, some of 
which have been suggested in the context of mean- field DLA In our earlier work on a variant of this model M , 

we employed the stratagem of cutting off the (x — x)P term when P{x) was under some threshold value P„. It is 
clear from the outset, however, that for this system such a resulting theory cannot give rise to a steadily-propagating 
solution, in contradiction with simulation. Furthermore, it is not clear how to justify this prescription (or any of the 
others that have been suggested) on a first-principles basis. Thus, we leave mean-field theory for the interim, and 
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study in detail some simpler, though nontrivial limits of the model, those of a very small number of individuals and 
the case of very small mutation rate \x. 

III. SIMPLE LIMITING CASES 

In this section, we discuss the simplest possible cases for the model as described above, that of populations of either 
one or two individuals or alternatively the case of extremely small mutation rate. As we shall see, this exercise allows 
us to determine exactly the nature of the long-time state and is furthermore quite useful in pointing out what one 
needs to do to make sense of the utter failure of simple mean field ideas. 

A. The Case of N = 1 

As an introduction to the case of 2 individuals, let us quickly see what happens for one individual. Although trivial 
to work out, it is actually a useful signpost on the road to the more interesting case of N = 2. 

For N = 1, we can write down the exact master equation for the Markov process. As there is only one individual, 
a configuration is labeled by the position of that one particle, so we can express everything in terms of f(n,t), the 
probability that the individual is at site n at time t. As the only way the individual can move to a new spot is by 
giving birth, which occurs at a rate n and then dying in favor of his child ( a 50% probability) so that his child can 
survive and, with probability /i, step to the left or right. Thus, the exact master equation reads 

f(n) = -\^nf(n) + [(n + l)/(n + 1) + (n - l)/(n - 1)] (4) 

or, in continuum terms, 

j{x) = -^{xf{x))" = \m{2f' + xf") (5) 

This is very reminiscent of the mean-field equation, Eq. [| except for the extra factor of 1/2 in the diffusion constant, 
but, more importantly the total absence of the troublesome xP term! The reason is clear: with only one individual, 
there can be nobody out in front of the pack out-reproducing everyone. In fact, the birth process itself is completely 
irrelevant in changing the configuration, as one of the two individuals existing post-birth is immediately tagged for 
death, leaving the configuration unchanged. 

As we shall soon see, this absence of the xP term is a very general result, and, in fact, the equation of motion 
for long times is, in the general case, exactly the same as for the N = 1 case, up to changes in the values of the 
coefficients. 

B. N=2: The Small fi Limit 

The case of N — 2 is less trivial. A general configuration is specified by the positions of both individuals, so that 
it constitutes a 2-dimensional field. However, things simplify tremendously in the small \i limit. The point is that if 
\x were exactly 0, then the population would quickly collapse, with both individuals finding themselves at the same 
location (most likely that of the fitter of the two). For small /i, every so often one of the two begins to wander 
away, but it can't get far before the population collapses again. So, for small [i, we need consider only 2 types of 
configurations: l)the two individuals on the same site, and 2) the individuals on nearest-neighbor sites. Call the 
probability the individuals are on the same site x, f(x). The probability that they are on x and x + 1 we call g{x). 
The master eqn. for the process is then: 

f(x) = ~nxf{x) + ^(1 - M ) (xg{x - 1) + xg{xj) + |((z + l)g(x) + (x - l)g(x - 1)) 

g(x) = ^(xf(x) + (x + l)f(x + 1)) - i(l - + 1) + x)g{x) - + 1) + x)g(x) (6) 

Notice we have forbidden the transitions from "g" to the state with a space between the 2 individuals, to conserve 
probability. We can of course combine terms above, but it is written so that each term in the / equation corresponds 
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directly to the same term in the g eqn. As we have argued, g is small, of order fi, so since we are working only to 
linear order in fi, we shall drop all terms involving ng. Then the above simplifies nicely to 



/(*) = - 



-fj,xf(x) + ^ (xg(x ~ 1) + xg(x)j 



)(x) = ^n(xf(x) + {x + l)f(x + 1)) - i ((x + l)g(x) + xg(x) 



(7) 



Notice now that there is an exact time-independent and ^-independent solution to these eqns., namely g = 2fxf = 
const. Of course, this solution is not normalizable, but as we will see, this is what the system is approaching for long 
times. This motivates us to guess that 



g(x) « 2uf(x) + fif'(x) (a x + nf'{x)B 



(8) 



or, equivalently (to the order we will need, (see below)) 



->,,/(,-! ^ ,,(,-) + ,/(,-) (Ci + ^) - </V)2> 



(9) 



Notice that this is an expansion in derivatives, which are assumed to decrease in magnitude with increasing order. 
The justification for this is the statement above, that / and g are approaching constants at long time, so the expansion 
is essentially one in inverse powers of t. We will see this explicitly later. It is also a large x expansion. In order 
to compute the velocity, we will need to keep terms up to x^ 1 in /', g' and up to x° in /", g" . We can relate the 
coefficients A and B to C and D by substituting the first expression in the second and expanding. We find A\ = — 2Ci, 
A 2 = -2C 2 , and B = -2C\ - 2D. 

The idea now is to substitute the above relation in the equation of motion for / and g. We note that / is proportional 
to but g has terms zeroth order in a. For this to be consistent with Eq. (||) relating g to /, we must make these 
zeroth order terms vanish. This gives us equations for the coefficients C and D. In more detail, 



g{x) = \ 
l 

~ 3 



(2x + l)g' (Ci + ^) + g"D + (x + l)(g> + C l9 ") + \g" 
xg\2C 1 + 1) + g'(d + 2C 2 + 1) + xg" '(2D + C x + h 



(10) 



Thus, C\ = —i, C 2 = — t, and D = 0, so that A\ — 1, A 2 — \, and B = |. Substituting the relation for g in term of 
/ in the eqn. for / yields 



f [Ai + —\ + f'B 



2A 2 



f [2Ai + — - 2 + f"(2B -A 1 + l) 



m = — 

fix 

~ T 

= 3 [f+xf] 

The last eqn. directly implies that the velocity v is since 



-^[2f' + A 1 f"] + if2f 



v = j f (dxxf{x) = ^ fxf + xif" = ^ 



(11) 



(12) 



where the last result comes from integration by parts and the fact that / is normalized to unity at lowest order. 
Furthermore, there exists a similarity solution for this / equation, (about which we will discuss more later) with x 
scaling like t, so that a long-time expansion is also a large-x expansion as advertised. Also, as advertised, the form 
of the equation of motion for / is exactly as in the N = 1 case. What we have therefore, is essentially a bound-state 
of the two individuals propagating together at constant velocity (on average) with the center of mass performing 
(inhomogeneous) diffusion. The squared-width of this bound-state, the expectation value of the squared-distance 
between the two individuals, is g/(f + g), which to the order we are calculating is simply 2u. 
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C. General N: The Small n Limit 



We can extend the calculation in the previous section to arbitrary N. To lowest order in jj,, we have to consider 
only those states where the individuals are distributed between 2 nearest-neighbor sites. We label the probability of 
being in the state with N — k individuals at x and k individuals at x ± 1 by (x), (0 < k < N/2). The probability 
of the state with all N individuals at x is f(x). If N is even, we also have the symmetric state gN/2{ x ), but, for our 
purposes it is enough to do the odd N case. Again, we can write down the master equation for the process, ignoring 
transitions out to other states, which are down by additional powers of \x. The equations of motion then read 



/(*) 

gt-i(x) 



l 



N + l 
1 

N + l 
1 

N+l 



N + l 



- N 2 »xf{x) + (N- l)x(g+{x) + g7(x)) 



l -N 2 ^xf{x) -(N- l)(2x ± l)gf(x) + 2(N - 2)xgf(x) 



— (N — k)k(2x ± l)g±(x) + (N - k - l)(n + l)xg± +1 (x) 
+ (N-k + l)(n-l)(x±l)g±_ 1 (x) 
N - 1 N + l 



(13) 



+ - 



2 2 
N-3 N + 3 



(2x±l)g^_ 1 {x) + 



(x±l)gj_ 3 (x) 



N — 1 N + l 



xg%_ 1 {x± 1) 



2 2 

where the middle equation holds for k = 2, . . . , (N — 3)/2. All the g's are again slaved to /, so we write 



9t{x) = rf{x) (a% ± ^ + + nf{x) (±Bl + ^) + M r 



(i)C fc 



(14) 



We have used the reflection symmetry to determine the ± signs in this equation. Essentially, they came from the fact 
that the ±l's in the equations of motion are down by 1/x from the leading terms. The procedure is exactly the same 
as before. The coefficients A\, B l k , and Ck are determined by the requirement that the g^r must all vanish to leading 
order. We find, after a mess of algebra, that 



A\ 



Bl 



Ck = 



N 2 



2k(N- 


k) 


N(N- 


- 2k) 


2k(N 


-k) 


N(k- 


1) 


2k 




N 





(15) 



2{N -k) 
N{N — k — 2) 
A{N - k) 
N-k-1 



N-k 



We see from this that the profile is asymmetric, with the left (trailing) side larger. However this asymmetry is 
proportional to 1 jx and so slowly decays in time. 

Now we plug all this into the equation of motion for /, finding that 



where the velocity, v, is given by 



/ = -vf + D(xff 
N(N - 1) 

v = M "2(ivTTy 



(16) 



(17) 



G 



and so increases (roughly) linearly with N, in accord with our general expectations above. The diffusion constant of 
the center of mass, D = 2 (n+i) 1S rou ghly constant with N. The (unnormalized) long-time similarity solution for / is 

/(*.«>- < 18 > 

We see that the spreading of the center of mass is linear in time. 

Unfortunately, we see that the g's are actually of order 0(fj,N), so the expansion is not uniformly valid in N. 
Nevertheless, it is clear that the form of the equation of motion is true to all orders in \i, with just the coefficients 
v and D varying. The task therefore is to compute these coefficients for arbitrary fi. For the case N — 2, we can 
generalize our analysis above to accomplish this, so we return to that case. 



D. Two Individuals, Arbitrary fi 

For the N = 2 case, we can straightforwardly generalize our small \i solution given above to the case of arbitrary 
fi. The relevant variables are f n (x,t), (n = 0,1,.. .), the probability of an individual at x and one at x + n. The 
equations of motion are: 



fo(x) 
h{x) 



fn(x) = 



-^xfo(x) + x(l - /x) lfn(x) + f n (x - n)} + + l)h(x + 1) + (x - l)h(x - 1)) 

2fi(xfo(x) + (x + l)f (x + 1)) - (1 + fJ,)(2x + l)h(x) 

+ [ X f"( X ) + ( x + l )fn{ x + 1) + Xf n (x - n) + (X + l)f n (x - 71 + 1) 

+ l»((x + 2)f 2 (x) + (x - l)f 2 (x - 1)) 
-(1 + n){2x + n)f n {x) + + n + l)f n +i(x) + (x - l)f n+1 (x - 1)) 

+ -n({x + n - l)f n -i(x) + (x + l)f n -!(x - 1)) 



(19) 



where in the last equation, n>2. Again, there is an exact steady-state solution, with all the /„ independent of x. If 
we choose the ansatz /„ = Ar n f 0} the equation for /„ reads 



= ~(2x + njAfor"- 1 (2(1 + n)r + fir 2 + m) 
which has the solution, (defining \ = l//i) 



X + 1 - vV + 2x 



(20) 



(21) 



For small /i, r ~ fi/2 <C 1, and r increases monotonically with fi. This is in accord with our expectation for a rapid 
falloff of /„ with n for small fi. For the maximum physical [i = 1, r = 2 — \/3 ~ 0.27. Plugging this form for /„ into 
either of the two remaining equations yields A = 4(1 — r)/(l — 3r). The normalization condition that J2n /« = 1 then 
determines /o = (1 — 3r)/(l + r). Note that fo and A are positive, as they should be, for all physical values of fi. We 
can then calculate the exact squared-width of the bound state, 



(22) 



Note that the first-order result derived above in the small \i limit is in fact exact to all orders! 

We can again generalize to the long-time limit where the derivatives of / are small, but non-zero. We write 



/„ = r n (Af + + C n f'Jx + D n fS) 



(23) 
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If we assume the expected form for /q: 



f = -vft + D(xf )"; (24) 

we can obtain a set of equations for the coefficients B n , C n and D n along with v and D by matching the coefficients 
of /o and its derivatives. After much algebra, we find 



B n = nA/2 



1 2 , (3 X -l)r_ , r 2 (2 x -3) 



C n = -T, n + ^ TT n 



8 4(1 -r 2 ) 2(1 -r) 2 (l +r)(l - 3r) 
2 ( X -3)r r 2 (2 X -l) \ 

D "=U n + 47I^ n+ 2(l-r) 2 (l + r)(l-3r)J A (25) 
v = /-i/3 
D = fi/3 

Just as was the case with the width, the coefficient of the drift v and diffusion D are given exactly by the linear in /i 
expressions found above! Such simple results cry out for a simpler derivation, which we now present. 



IV. MOMENT EQUATIONS 

As we have seen, the dynamics of our discrete evolution model gives rise to a pulse of a typical width; the center of 
this pulse propagates at long times with some velocity and also diffuses. Unfortunately, the methodology used to obtain 
these results (essentially studying the full master equation) cannot be used to obtain the macroscopic parameters of 
the velocity and diffusion constant in general, and also do not indicate anything about the pre-asymptotic state. What 
we will see now is that one can re-derive the same picture for the long-time limit of the model at arbitrary parameters. 
This derivation proceeds by working with moment equations for the probabilistic process. These moment equations 
are valid at all times, and so also shed light on the short-time dynamics of the system. 



A. The Width of the Pulse 



Staying for the moment with the N = 2 case, let us focus on the average squared- width, < n 2 >. We wish to 
calculate < n 2 >, the rate of change of the squared- width with time. Given the definition of < n 2 > in terms of the 
fn's, (see Equation (p2|)) we can obtain what we want from the equations of motion of the / n 's. We find 

2 , 4 , 

< n 2 > = --^in 2 /„(i) + ^^x.fn(x) 

n n 

2 4 

= — - < in > +-/_* < x > (26) 
o o 

where x = x+n/2 is the average position of the two individuals making up the configuration f n (x). If we may factorize 
the expectation value < xn 2 >«< x >< n 2 >, (a step we can justify for long times from our explicit solution above) 
we get 

2 4 

< n 2 > = -- < x X n 2 > < x > (27) 

o o 

which immediately yields the steady-state value of lu = 2/i, which we obtained in the previous section. 

This argument can be easily extended to arbitrary N. We now define the squared-width, uj, as the mean squared 
distance between pairs of individuals, taken over all N(N — l)/2 pairs. 

v > i>j 

and so its ensemble average is 
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< OJ > = 



N(N - 1) 



> C i>j 



(29) 



where the first sum is over all configurations C, weighted by their probability Pc- We are thus considering a 2-point 
function, instead of the full N-individual joint probability distribution. We can calculate to by considering the effect 
of each possible event on ui. We find, summing over all events where i is born and j dies, 



< u> > = 



N(N~ 1)(N + 1) 



(I -/,) ( £ [(x k 




Xi) 2 - (xk - Xj) 2 ] + - (xi - Xjf 



l) 2 - (x k - Xjf] + l~{Xi 



+ l) 2 - (x k - Xj) 2 ] +l-(xi- Xjf 



After some algebra, we get 



< lu > = 



2N 
N+l 



< xoj > 
N 



+ < C 3 > +n < x > 



(30) 



(31) 



where x is the average x value of the population, (l/AOX^ 2 ^ an d C3 is the asymmetry, (l/-^)X)i( x i — ^) 3 - This 
asymmetry automatically vanishes in the case of 2 particles, so the equation reduces to the Equation (27) in that 
case. We have already seen that the average asymmetry vanishes at long times like 1/t, (an argument which easily 
extends to all orders in /1), so that upon factorizing the expectation value < xoj > as above, we get that the width 
asymptotically approaches /J.N. 

We can understand simply the origin of each of the three terms in the equation of motion for < uj >. The fi term 
comes from the fact that mutation acts to increase the width. Since only babies can mutate, the effective mutation 
rate is proportional to fix. The asymmetry term expresses the effect of selective pressure. If there are particles that 
lag behind, giving rise to a negative asymmetry, they are quickly killed off, decreasing the width. The last term, 
< xoj >, is more subtle. It expresses the fact that the population tends to collapse even in the absence of selective 
pressure. Death always acts to reduce population variability, whereas birth does not increase it. This effect is crucial 
for the stabilization of the width against the effects of mutation at long times. The rate associated with this effect 
is 2x/(N + 1). As the velocity scales linearly in N, however, (which we have proved to lowest order in /1 above and 
which we demonstrate below to all orders) the time scale for the asymptotic width dynamics is essentially independent 
of N. Unfortunately, this does not mean that the asymptotic state is necessarily approached in a time which is N 
independent; we will see later that this issue is actually considerably more complex. 



B. The Drift Velocity 

It turns out that not only can we derive the width to all orders in /1, we can also simply relate this width to the 
velocity of the average mean position, x. 
The average mean position is given by 

<x>=^Y, p cll xi ( 32 ) 

C 1 

Now a little reflection suffices to convince oneself that mutation does not contribute to the instantaneous velocity, 
since the probabilities for mutating left and right are identical. Thus the entire velocity is due to the collapse events. 
Then, as above, summing over all events where i has a baby, resulting in the death of j, 

<*> = N J +1) Y l 12X' PcXi{Xi ~ Xj) 
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c 

N - 1 



2(JV + 1) 



< uj > (33) 



Thus, the velocity is just exact proportional to < uj >, the average squared- width! This is physically reasonable: the 
velocity arises from the collapse of the population onto the more fit members. The wider the distribution, the faster 
the average fitness increases. Now, we found above the asymptotic value of < uj > to be /iiV, yielding the striking 
result that the velocity is exactly linear in /i, and given by our lowest order result, Eq. ( |l7| ) found above! 

Turning back to the simulation results in Fig. 2, we see that the results are indeed consistent with a velocity 
converging to the theoretical value, Eq. (fL7|). Extending the run to longer times confirms this result. Similarly, we 
see the squared-width tracks with the velocity as predicted, increasing toward the expected value /iN. 

We can also calculate in a similar manner the diffusion constant D of the center of mass motion. From Eq. ( |l6| ) 
above, one can deduce that d/dt (< (x) 2 > — < x > 2 ) = 2D < x >. Calculating this quantity directly via the 
moments, we find 

d / , sO 9N 1 f N-l \ 

- (< {xf >-<x> 2 ) = j^-^ I < C 3 > + — ^— <xoj>+[Kx >\ (34) 
Since, at long times, < C3 > decays to zero, and < oj > approaches [iN , we get that 

2(N+1 



again reproducing our lowest order result exactly. 

So, we have been able to show that the asymptotic state of the evolution process is one in which the width saturates 
to a fixed value, leading to a pulse propagating at constant velocity in fitness space. The exact values of the velocity 
and width have been determined. What we still need to discuss is the failure of mean-field theory and the related 
issue of the transient behavior before the asymptotic state sets in. We will return to these issues shortly. Now, we 
briefly detour to comment on the flat landscape case. 



V. THE FLAT LANDSCAPE AND OTHER VARIANTS 



All the analysis above works as well (and more simply) in the case of the flat landscape; i.e., where the birth rate is 
independent of x. A variant of this flat-landscape model in which not only babies mutate, but every individual does, 
was considered first by Zhang, et. al || and exactly solved by Meyer, Havlin, and Bunde fli]|| . In this variant, the 
relevant parameter is the ratio of the mutation rate to the birth rate, which we may take to be unity. This ratio we 
label fj,, to distinguish it from /x, the probability that a baby mutates. The context in which this variant was analyzed 
was a problem of population migration, so that the "mutation" occurs in physical space, not "fitness" space. One can 
consider also an intermediate case, where everyone mutates (or migrates) and the birth rate depends on position, (the 
"emigration" problem, where people look for "greener" pastures). All these variants have many features in common. 
As we shall see, the major differences between the models are l)the nature of the "drift", or velocity; and 2)the 
time-scale for the wave-function to reach its asymptotic width. 

The flat landscape model of Zhang, et. al. and Meyer, et. al. is sufficiently simple to allow for a complete analytic 
solution, which Meyer et. al. present in a lovely piece of work One essential difference in the flat landscape case 
is that there is no drift. The equation of motion for the center of mass in the long-time limit is exactly the diffusion 
equation: 

./ = Df" (36) 

where the diffusion constant D — fl/2 is exactly that of a single individual. The explanation for this fact, already 
noted by Zhang, et. al., is that the collapse event leaves the parent unaffected, so that there is at long times one 
"ur-individual" who just propagates as if nothing ever happened. The localization of the wavefunction then implies 
that all the diffusion of the center of the mass is due just to the diffusion of this "ur-individual" . Comparing this to 
our result we see that this parallels our results for the diffusion constant in the "tilted" case. There we found that 
the diffusion constant was N[xx/2(N + 1). This is just the effective diffusion constant for a single individual, since in 
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the babies-only model, a diffusion step only occurs if the baby is not immediately killed, so that the effective rate of 
mutation is fiN/(N + 1) times the birth rate x, giving rise to a diffusion constant as noted above. The fact that only 
babies mutate is not an essential change since the dynamics would be the same if only the parent mutated, so again 
only the motion of the "ur-individual" need be considered. 
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FIG. 3. a) Mean birth-rate for population of N — 80 individuals as a function of time in the "emigration" model, with 
mutation rate /i = 5., averaged over 100 realizations. b) Width-squared from same simulation. 

Of course, the drift must vanish in the flat-landscape case, by virtue of the left-right symmetry In the intermediate 
model where all individuals can mutate, but the landscape is tilted, it is straightforward to see that the velocity is 
proportional to 1/x, the current average position, so that the position should scale as \ft. This is related to the fact 
that the width of the wave function decreases as 1/x. The point is that a constant mutation rate is not effective enough 
to combat the ever-increasing suppression of variability due to the ever-increasing birth/death rate. This behavior is 
apparent in the simulation data of this model for N — 80 presented in Fig. 3. We will return to a discussion of this 
"emigration" model in Sec. VII. 

As already mentioned, the third order moment C3 in the width equation represents the effect of selective pressure. 
Is is therefore absent in the flat landscape case; instead, the moment equation for the width squared, u), reads 

<c> = 2(-^+m) (37) 

Also missing, naturally, are the factors of x. Thus the moment equation closes, allowing for an exact solution, identical 
to that found by Meyer, et. al. 

The last major difference worth noting is the time scale for the asymptotic width dynamics. We saw how in the 
tilted case the time scale is x/N. The factor x is the effective birth rate. In the flat landscape case, the birth rate 
is constant, so the time scale is fixed and or order l/N. Thus stabilization of the width in the flat landscape case 
always takes place after O(N) time. 



VI. MOMENT EQUATIONS AND MEAN-FIELD THEORY 

The equations of motion we have derived for the width and velocity not only allow us to calculate the asymptotic 
behavior of these quantities. They also allow us to make contact with the mean-field theory we presented in section 
IIB. The point is that one can derive parallel equations for these moments directly from the mean- field (MF) theory. 
By comparing the exact equations with the approximate ones obtained from the mean-field theory, we can obtain 
some insight into the breakdown of the mean-field theory, and how to improve upon it. 
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The MF equation for the velocity is easy to obtain. The mean of x is 



5>P(n) (38) 

Using the MF equation of motion for P, Eq. (^), we obtain 

< x >MF = GfF (39) 

where C^ F is the second cumulant, < n 2 >mf — < n >L. This second cumulant is, to leading order in 1/N, 
just lo/2. Thus, comparing with the exact equation, Eq. (|33|), we see that the MF result is exact up to trivial 1/N 
corrections. 

We can obtain in a similar manner an equation of motion for the second cumulant. We find 

Cf F = Cf F + M < x > MF (40) 

Comparing with the exact equation for u>, Eq. (^lj) we see the equations correspond except for the absence of the 
stabilizing term — < XC2 > /N. Now of course this term is formally of order 1/iV, and so the MF cannot be expected 
to reproduce it Nevertheless, we have seen that this is the term responsible for the saturation of the width at 
long times. The point is that < C2 > is O(N) at long times and so the term is in truth not negligible. 

The question at this point is, how long is long? When does this and other formally small terms become relevant? 
The answer is very quickly, on a time scale of 0(1). The reason for this is the blowup of MF we alluded to in section 
2. We can see this blowup directly by solving exactly the MF equation of motion, Eq. To do this, we define a 
generating function for the moments, Z(a,t), 

Z(M)^]TP(M)e a " (41) 

n 

It is actually more convenient to work in terms of the function F = In Z, which is a generating function for the 
cumulants. Plugging the definition of F in the equation of motion for P, we find 



F(a,t) = (1 + ^cosh(a) - fi)-^-F(a,t) - ^f- 

da da 

We can solve this equation by means of defining a new variable s such that 

ds 1 



a=0 



da 1 + [i cosh(a) — \i 



(42) 



(43) 



so that 



V ; VI - 2/i V 1 - V 1 - 2/itanh(a/2) J y ' 

for /i < 1/2 and 

s(a) = 2 arctan - ltanh(a/2) S ) (45) 

for /i > 1/2. In terms of s, the equation for F reads 



■ , , d , , dF 
Hs,t) = -F M -- 



(46) 



s=0 



whose solution is 



F(s,t)=F (s + t)-F {t) (47) 

where Po(s) = Fo(a(s)) is the generating function for the cumulants of P(n,t = 0). For example, if we start with a 
zero width pulse located at no at time t = 0, Fq(s) = noa(s). Then, translating back to our original variable a, we 
find (for our zero width pulse) 
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F{a,t) = n (s- 1 (s(a) +t) - s" 1 ^)) 



(48) 



The key feature of this solution is that = s(a — oo) is finite, so that the inverse function s 1 is singular at Soq. 
This causes F to diverge at time t* — Sqo. Explicitly, 

1 Tsf^T arctan(V2/x - 1) /x > \ 

so that i* diverges as — In// for /x <C 1 and decreases monotonically to tt/2 for the maximal fx = 1. For example, the 
first moment of P, the mean position, is given by 



dF 

< X >MF= -J— 

da 



(1 lj-tanh 2 (^))cosh 2 (^) ^ ^ 2 



' ' i (50) 



(1- i tan 2 (¥))co S 2(M) M > 



where x = vP 2/i|. The divergence of this quantity at i = t* is apparent. 

In fact, all nonzero cumulants of P diverge at t* , with the jth cumulant diverging as (t* — t)~ J . This makes an 
analysis of the precise manner of the breakdown of MF rather subtle. Consider the equation of motion for lo, Eq. 
(pl|). It is true that < x > ui/N is diverging as t — ► t* as (f* — i) -3 , but then again the C3 term in the equation is also 
diverging in exactly the same manner. As it is clear that MF must break down before t = t*, it must be higher-order 
moments that are responsible. These higher-order moments are responsible for driving the low-order moments to 
diverge, and if they are cut off somehow by the 1/N terms, the effect will eventually propagate down to the low-order 
moments. 

It is difficult at this stage to pinpoint precisely the workings of this mechanism however. At issue is whether there 
exists an intermediate time regime at very large N where neither MF nor the asymptotic regime with its constant 
velocity are accurate descriptions of the dynamics. To indicate the nature of the problem, we present in Fig. 4 a 
graph of the time-development of the squared- width, < lo >, as a function of time as measured from simulations for a 
number of different iV's. The initial time regime with its MF-like behavior is clearly evidenced for short times. This 
rapid rise in the velocity is then abruptly cut off. The velocity then slowly increases to its asymptotic value. The 
time-scale of this slow rise, and its dependence on N are not at this point clear to us. We hope to return to these 
issues in a future work. 



VII. THE "EMIGRATION" MODEL 



We now examine the breakdown of MF for the variant model where everyone mutates, the "emigration" model. As 
we indicated in Section VI, the asymptotic state is characterized by a velocity decreasing as 1/Vt, The short-time 
dynamics is controlled by the following MF equation: 

P(n) = (n - n)P(n) + -p, [P(n + 1) - 2P(n) + P(n - 1)] (51) 

We can also solve this MF exactly. The equation for F(a) now reads 

(52) 

For an initial population located at no, so that F(a, t = 0) = nga, one can verify that 

F(a, t) = jj,( sinh(a + t) — sinh(i) — sinh(a)^ + noa (53) 

satisfies the MF equation. From this solution, one can read off that the velocity v = fx sinh(i) grows exponentially 
in time. Thus, for this model, MF must break down by the time ln{N/p) since at this time the MF value of 
< x >< lo > /N is competitive with the C3 term as given by MF. However, as with the only-babies-mutate model 
analyzed in the previous section, the breakdown of MF appears to be more complex. Simulations indicate that in fact 
MF breaks down before this time and in particular when < x >< lo > /N is still small compared to the other terms. 
We present in Fig. 5, simulations of the width as a function of time for various TV. Again, for very short times, we 



F(a) = -F(a) 



/2(cosh(a) — 1) — — 



dF 
da 
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see a regime where MF is accurate. The growth of the width is then suddenly cut-off and the width starts decaying. 
The cutoff is seen to take place at times earlier than that indicated by the simple argument above, due to the role of 
higher-order moments, similar to what we found for the only-babies-mutate model as discussed above. 
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FIG. 4. Squared- width as a function of time for populations of N = 100 (circle) , N = 200 (diamond) , and N = 400 (triangle) 
individuals, with mutation rate fj, = .5, averaged over 160, 640 and 160 realizations respectively. 

In our previous work Q, we presented an analysis of a cut-off version of the MF for this everyone- mutates model. 
The cut-off amounted to killing the n — n term for those sites where P(n) was less than some threshold value [g). 
This model possessed an asymptotic state which propagated at constant velocity. In fact, starting from a distribution 
localized on a single site no, the velocity in the cut-off theory increases monotonically to this asymptotic value. We 
see that this cutoff theory seriously misrepresents the actual physics, even at large N. In essence, this form of cutoff 
is too weak. We have seen that the corrections to MF are in fact stabilizing in nature, whereas the simple truncation 
of this cutoff just makes the growth at the leading edge marginal. The question of whether a more physically correct 
cutoff can be constructed is intimately related to the question of the precise nature of the breakdown of naive MF 
and the rapidity of approach to the asymptotic state; these questions we have chosen to postpone to future work. 

VIII. CONCLUSIONS 

The purpose of this paper was to present a detailed analysis of a simple model for evolution on a smooth landscape. 
This model was originally motivated by experiments on RNA viruses which showed persistent increases in fitness 
without the "punctuation" that should be present if the fitness landscape had significant roughness thereby giving 
rise to local maxima. Our results indicate that the system does settle into a asymptotic state characterized by a 
constant rate of fitness growth. This state is dominated by a balance between mutation and the the loss of variability 
within the population due to death. This latter effect, although formally weak in a large population, is crucial in 
rescuing the population from the runaway increase in birth-rate seen in the mean- field limit. 

Our model exhibits certain generic features which should be sufficiently robust that it is reasonable to hope that 
they apply also to the biological system which motivated these investigations. The most striking is the role of the 
population size on the dynamics. Since the iV-independent mean-field treatment breaks down so spectacularly after 
a very short time, we may expect that the population size is a very relevant parameter. The other prediction of our 
model that bears checking in the experimental system is the behavior of the population variance. We have seen that 
if the initial variance is large, there is a rapid initial collapse of variance. The variance then asymptotes to a fixed 
quantity and does not exhibit diffusional broadening in time. Also, the fact that the rate of improvement of fitness 
vanishes for small N is very relevant for modelling the process of "genetic bottlenecking" jl3| used to create the initial 
population in the experiment. 

Our investigations have implications beyond the narrow confines of the model analyzed herein. The breakdown of 
mean-field theory exhibited by our model is shared by diffusion-limited aggregation. There also it is clear that this 
breakdown is the result of the neglect of correlations, in this case correlations between the walker density and the 
density of the aggregrate. This correlation is most apparent for particles who manage to outrun the typical extent of 
the aggregate. Whereas the ensemble-averaged walker density is quite large at the position of this outlying particle, 
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the actual walker density is very small. The neglect of this correlation leads to a blowup of the mean-field exactly 
parallel to the blowup exhibited by the mean-field studied herein. The failure of simple cutoffs to correctly capture 
the asymptotic dynamics leads to the concern that similar cutoffs employed in the study of mean-field DLA might 
also misrepresent the dynamics of the system. 

Our methods might also prove useful in the study of the behavior of genetic algorithms Q, where presents a 
similar system of evolving populations. Questions such as the speed of approach to the optimal state, the dependence 
on population size and the behavior of the population variance lie at the heart of understanding the working of these 
algorithms. 
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